home *** CD-ROM | disk | FTP | other *** search
/ Pascal Super Library / Pascal Super Library (CW International)(1997).bin / MATH / MATH1 / TSTBES.PAS < prev    next >
Pascal/Delphi Source File  |  1985-04-03  |  2KB  |  83 lines
  1. program tstbes;        { -> 344 }
  2. { test the bessel function }
  3. { the Gamma function is included }
  4.  
  5. var    done    :boolean;
  6.     x,ordr    : real;
  7.  
  8.  
  9. function gamma(x: real): real;
  10. const    pi    = 3.1415926;
  11.  
  12. var    i,j    : integer;
  13.     y,gam    : real;
  14.  
  15. begin        { gamma function }
  16.   if x>=0.0 then
  17.     begin
  18.       y:=x+2.0;
  19.       gam:=sqrt(2*pi/y)*exp(y*ln(y)+(1-1/(30*y*y))/(12*y)-y);
  20.       gamma:=gam/(x*(x+1))
  21.     end
  22.   else        { x<0 }
  23.     begin
  24.       j:=0;
  25.       y:=x;
  26.       repeat
  27.     j:=j+1;
  28.     y:=y+1.0
  29.       until y>0.0;
  30.       gam:=gamma(y);        { recursive call }
  31.       for i:=0 to j-1 do
  32.     gam:=gam/(x+1);
  33.       gamma:=gam
  34.     end    { x<0 }
  35. end;        { gamma function }
  36.  
  37. function bessj(x,n: real): real;
  38. { cylindrical Bessel function of the first kind }
  39. { the gamma function is required }
  40.  
  41. const    tol    = 1.0E-4;
  42.     pi    = 3.1415926;
  43.  
  44. var    i        : integer;
  45.     term,new_term,
  46.     sum,x2        : real;
  47.  
  48. begin    { bessj }
  49.   x2:=x*x;
  50.   if (x=0.0)and(N=1.0) then bessj:=0.0
  51.   else if x>15 then { asymptotic expansion }
  52.     bessj:=sqrt(2/(pi*x))*cos(x-pi/4-n*pi/2)
  53.   else
  54.     begin
  55.       if n=0.0 then sum:=1.0
  56.       else sum:=exp(n*ln(x/2))/gamma(n+1.0);
  57.       new_term:=sum;
  58.       i:=0;
  59.   repeat
  60.     i:=i+1;
  61.     term:=new_term;
  62.     new_term:=-term*x2*0.25/(i*(n+1));
  63.     sum:=sum+new_term
  64.   until abs(new_term)<=abs(sum*tol);
  65.   bessj:=sum
  66.  end    { if}
  67. end;    { bessj }
  68.  
  69. begin        { main }
  70.   done:=false;
  71.   repeat
  72.     write('Order: ');
  73.     readln(ordr);
  74.     if ordr<-25.0 then done:=true
  75.     else
  76.       begin
  77.     write('X: ');
  78.     readln(x);
  79.     writeln('J Bessel is ',bessj(x,ordr))
  80.      end
  81.   until done
  82. end.
  83.